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Abstract 


This report discusses the development and application of two alternative strategies in the form of 
global and sequential local response surface (RS) techniques for the solution of reliability-based 
optimization (RBO) problems. The problem of a thin- walled composite circular cylinder under 
axial buckling instability is used as a demonstrative example. In this case, the global technique 
uses a single second-order RS model to estimate the axial buckling load over the entire feasible 
design space (FDS) whereas the local technique uses multiple first-order RS models with each 
applied to a small subregion of FDS. Alternative methods for the calculation of unknown 
coefficients in each RS model are explored prior to the solution of the optimization problem. 
The example RBO problem is formulated as a function of 23 uncorrelated random variables that 
include material properties, thickness and orientation angle of each ply, cylinder diameter and 
length, as well as the applied load. The mean values of the 8 ply thicknesses are treated as 
independent design variables. While the coefficients of variation of all random variables are 
held fixed, the standard deviations of ply thicknesses can vary during the optimization process as 
a result of changes in the design variables. The structural reliability analysis is based on the first- 
order reliability method with reliability index treated as the design constraint. In addition to the 
probabilistic sensitivity analysis of reliability index, the results of the RBO problem are 
presented for different combinations of cylinder length and diameter and laminate ply patterns. 
The two strategies are found to produce similar results in terms of accuracy with the sequential 
local RS technique having a considerably better computational efficiency. 

1. Introduction 

The desire to increase structural efficiency while maximizing reliability and robustness has 
fueled the rapid growth of non-deterministic approaches for more than a decade, and has led to 
the development and application of many techniques under the general category of reliability- 
based optimization (RBO). Numerous techniques have been developed to assess the structural 
reliability of a component or a system, which can be broadly categorized into two groups: (1) 
random sampling methods, and (2) analytical methods. In all cases, the uncertainties associated 
with structural modeling, material properties, loading condition, etc. are captured using a variety 
of techniques that rely on probabilistic or fuzzy treatment of random variables. 

Prevalent among the random sampling methods are the traditional Monte Carlo simulation and 
its variants such as adaptive importance sampling (e.g., Wu 1994) and robust importance 
sampling (e.g., Tomg et al. 1996). These methods have been developed to reduce the simulation 
cycles by adjusting the sampling domain to the important region near the most probable point 
(MPP) of failure. 

Analytical methods such as the first- and second-order reliability methods (FORM and SORM) 
are based on the formulation of the limit state function that defines the surface separating the 
failure and safe regions of the design space. These techniques are used to determine the 
reliability (safety) index represented by the shortest distance from the origin of the reduced 
coordinate system to the limit state surface. A commonly used procedure for the calculation of 
safety index is that developed by Hasofer and Lind (1974) and expanded by Rackwitz and 
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Fiessler (1978) to include the effect of distribution type as well as the mean value and standard 
deviation of each random variable. 

The calculations of reliability index as well as the probability of failure often require repeated 
evaluations of the limit state function as a result of perturbations in the random variables in 
search of MPP. This repeated evaluation of the limit state function and its sensitivity derivatives 
are highly dependent on the complexity and computational intensity of the underlying failure 
analysis as well as the number and distribution type of random variables considered. In order to 
reduce this computational burden, it is possible to develop an accurate analytical model as a 
surrogate for the more exact failure analysis. This approach is made possible with the help of 
response surface methodology (RSM). 

Bucher and Bourgund (1990) used an adaptive interpolation scheme to represent the system 
behavior by a response surface (RS) model, which was then combined with importance sampling 
to obtain the desired reliability estimates. Liu and Moses (1994) combined the sequential RSM 
with Monte Carlo importance sampling to develop a reliability analysis program for aircraft 
structural systems; whereas, Liaw and DeVries (2001) developed a reliability-based optimal 
design process by integrating reliability and variability analyses with optimization design 
processes using RSM. 

The research presented in this report examines the use of RSM in reliability-based optimization 
of composite structures through the development of two alternative strategies, one dubbed the 
global response surface technique (GRST) and the other sequential local response surface 
technique (SLRST). 

These two strategies are applied in RBO of a thin-walled composite circular cylinder under axial 
compression with the results used to compare the accuracy and efficiency of these techniques. 

The remainder of this report is organized as follows. In Section 2 a brief overview of structural 
reliability is presented. This is followed by a detailed discussion of the RBO example problem 
in Section 3, which includes the description of GRST and SLRST and corresponding results. 
Reliability sensitivity analysis results are presented in Section 4 followed by the concluding 
remarks in Section 5. The description of the shell buckling analysis is given in Appendix A 
while that for the selection of design points for SLRST implementation is given in Appendix B. 

2. Overview of structural reliability 

In its fundamental form, the probability of failure of a structural component may be expressed as 
P f -P(g(X)s 0) (1) 


where X - {Xj, X2,—,X n \ T is an n-dimensional vector of random variables an dg(X) is the limit 
state function describing the failure criterion, such that g < 0 represents failure, g > 0 represents 
safety, and g = 0 represents the limit state surface separating the safe and failed regions. The 
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failure probability is found by evaluating the multiple integral of the joint probability density 
function fx(x) over the failure region, Q expressed as 

Pf B Sf -fafx( x l’X2>- ,x n )dx l dx 1 ...dx n (2) 


For problems involving multiple random variables, the integration of f x (x) is in general very 
difficult. Hence, the probability of failure is usually estimated using variety of approximation 
techniques including those commonly known as the random sampling methods (e.g., Monte 
Carlo simulation), analytical methods (e.g., FORM), or hybrid methods (e.g., AMV+ combined 
with AIS, Wu 1994). 


In applying the standard time-invariant FORM, the limit state g(X ) is transformed to standard 
normal space g(u ) with MPP found by solving the optimization problem 



S. T. g(u) = 0 


( 3 ) 


where the distance to MPP is defined as the reliability or safety index fi = |« |, which is then 
used to estimate the probability of failure (for /3 > 0) as 

P f - $(-/3) (4) 

where d> is the cumulative distribution function of a standard normal random variate 


The optimization problem in Eq. (3) may be solved using a variety of mathematical 
programming techniques such as that developed by Shinozuka (1983) or other more tailored 
techniques including those developed by Hasofer and Lind (1974) and extended by Rackwitz and 
Fiessler (1978) to problems involving non-normal random variables. In this research, we used 
the standard time-invariant FORM to find /?. 

One advantage of analytical over simulation-based techniques is that they facilitate the 
calculation of probabilistic sensitivity factors represented by the partial derivatives of (5 with 
respect to the mean and standard deviation of random variable X, as (Madson et al. 1986) 



(5-a) 


(5-b) 
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While the sensitivity factors in Eq. (5-a) can help measure the influence of each random variable 
on /?, those in Eq. (5-b) quantify the effect of parametric uncertainty on component reliability. 

3. Reliability-based structural optimization problem 

The goal of RBO is to achieve structural efficiency while fulfilling the requirement for safety. In 
this regard, parameters with uncertainty are treated as probabilistic random variables with each 
defined by its mean value, standard deviation, and probability distribution. Additionally, each 
failure constraint is formulated in terms of probability of failure or the corresponding safety 
index. 

The design problem considered here is a thin-walled circular cylinder under uniform axial 
compression as shown in Fig. 1, with the loaded edges clamped. 
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Fig. 1 Structural model of thin- walled composite cylinder and loading condition 

The cylinder wall consists of a 16-layer symmetric laminate with each layer made of 
unidirectional carbon/epoxy tape with statistical properties shown in Table 1 . In addition to the 
material properties, the thickness and orientation angle of each layer on one side of the plane of 
symmetry are also treated as independent random variables with the mean thickness in layers 1 
through 8 treated as design variables. The cylinder length (L) and inner diameter (D) are allowed 
to have three specific mean values as indicated in Table 1. Completing the list of random 
variables is the resultant axial load, P a . In this problem, all random variables have normal 
distribution with assumed coefficients of variation (COV) given in Table 1. 

The RBO problem in this example is formulated as 

Min £[W(X)] 

S.T. h * /S min (6) 

y! zY,z Y“, 


and 


i - 1,2,. ..,8 


where the objective function represents the expected or mean cylinder weight and represents 
the reliability index associated with axial buckling load, which is assumed to be the only mode of 
failure with the method of analysis described in Appendix A. The target component reliability 
index is denoted by /3 min . The design variables T, through Y % correspond to a subset of random 
variable vector, X and are constrained by their respective lower and upper bounds. Although the 
COV of each ply thickness is held fixed, the corresponding standard deviation is subject to 
change during the optimization process as a result of changes in mean ply thicknesses (i.e., 
design variables). 


Table 1 Statistical properties of random variables 


Random Variable 
No. Definition Unit 

Mean 

Value 

COV 

(%) 

1 

Ei 

psi 

18.0e6 

1.33 

2 

E 2 

psi 

1.35e6 

1.67 

3 

V I2 

— 

0.226 

2.0 

4 

^12 

psi 

0.543e6 

2.0 

5-12 

^piy 

in 

[y,/y 2 /.../y 8 ] 

5.0 

13-20 

0ply 

deg 

[d l ie 2 /.Jd s ] 

00 * 

21 

D 

in 

10 s E[D] =s 20 

5.0 

22 

L 

in 

10 <; E[L\ 20 

5.0 

23 

£ 

lb 

0.10e6 

5.0 


“Standard Deviation 


The optimization problem in Eq. (6) is solved using the method of sequential quadratic 
programming in the optimization software, DOT (1999) based on two alternative strategies for 
evaluation of p b . As will be described in the following section, these strategies take advantage of 
response surface methodology to increase the efficiency of component reliability evaluation, 
which is the most time-consuming part of reliability-based design optimization process for the 
types of problem such as the one considered here. 

3.1 Global response surface technique (GRST) 

In applying GRST, a single second-order RS model in n factors expressed as 

P cr (X) -ao + latXi+i 2 qjXjXj (7) 

i-1 i- 1 7-1 


and denoted as the global response surface model (GRSM) is used to estimate the axial buckling 
force, P cr over the entire feasible design space. This is done in lieu of calculating the "exact" 
P cr from a more costly buckling analysis described in Appendix A. Equation (7) is a fully 
quadratic model with (n + l)(n + 2) / 2 coefficients. With the cylinder axial buckling affected by 
random variables 1 through 22 in Table 1, Eq. (7) will be a function of 22 factors for a total of 
276 unknown coefficients. 
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To evaluate the buckling reliability index, using FORM, we approximate the limit state 
function as 

g(X) - g(X) = P cr - P a (8) 

where P cr is calculated from Eq. (7) and P a , representing the applied axial force, is defined as 
random variable X 2} . Since g(X) is a nonlinear function of random variables, the evaluation of 
P b becomes an iterative process within the RBO Solution block in the flowchart of optimization 
procedure shown in Fig. 2. 



Fig. 2 Flowchart of optimization procedure with GRST 

Thus, the main challenge in implementation of GRST is the development of an accurate GRSM 
for P cr . For generating the necessary response population in fitting a second-order RS model 
commonly used methods include experimental design (e.g., central composite design, Box- 
Behnken) and random sampling (e.g., Latin Hypercube). In this problem we used the method of 
least squares to estimate the unknown coefficients in Eq. (7) based on a population of response 
samples generated using Monte Carlo simulation (MCS), which is more robust but somewhat 
less efficient than the stratified random sampling and some of the available experimental design 
methods. 

In performing MCS, we assumed each random variable has a uniform distribution with the 
corresponding lower and upper limits found as 

X], x}-(l± Sx .)v X j, (9) 

where gx j represents the bound increment (BI) on ith random variable with all values shown in 

Table 2. The values of BI for material properties extend the range of GRSM to a minimum 
of ju± 3a for each property. The large BI for ply thickness set the min and max values to 0.0025 
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in. and 0.0075 in. For d p i y * 0 , the selected BI extends the domain of GRSM by considering 

variability and uncertainty in ply orientation angle. Finally, the large Bis for length and diameter 
allow the use of GRSM for the range of cylinders identified in Table 1. 

Table 2 Bound increments on random variables 


Random Variable 

SXj ( % ) 

E\ 

4 

E 2 

5 

Vw 

6 

G, 2 

6 

^ply 

50 

fyly 

26 

D 

50 

L 

50 


A deterministic computer program based on the procedure described in Appendix A is used to 
calculate the "exact" buckling load with the 12th degree Legendre polynomial as the 
interpolation function. Each analysis run takes approximately one minute and fifteen seconds of 
CPU time on a 900MHz Sun Fire 480R. For the subsequent numerical analyses, we considered a 
cylinder with mean ply angles and laminate stacking sequence specified as (±45 / 9 O 4 / +45 ) . 

For response surface analysis, we defined Eq. (7) as a multiple linear regression model with the 
least-squares approximations of the unknown coefficients found using a combination of 
GLMMOD and REG procedures in SAS version 8.0. 

Several candidate models based on increasing number of response samples were developed and 
compared prior to selecting the most appropriate GRSM for reliability-based optimization. The 
results of statistical analysis are summarized and shown in Table 3 with the last two rows in the 
table showing the results for a linear model (i.e., 3000L) and that excluding all the interaction 
terms (i.e., 3000LQ). 

In terms of RMSE adjusted R 2 , and COV, the model based on 500 samples appears to have a 
slight advantage over the rest. However, according to Freund and Littell (2000), the presence of 
significant round-off errors could cause the predicted residual sum of squares (PRESS) to be 
drastically different from the sum of squares of actually computed residuals (SS Error). In this 
case, this ratio is very large for the 300-sample case but quickly drops as the sample size is 
increased. The results in the last two rows indicate that in this problem, the interaction terms are 
very important and their exclusion would greatly compromise the accuracy of the resulting RS 
model as indicated by the much larger RMSE value. 

In addition to the statistics shown in Table 3, we also used the variance inflation factor (VIF) to 
detect the presence of multicollinearities among the 22 independent variables. In this case, the 
VIF values were found to be in the range of 1.003 to 1.012 indicating no multicollinearities 
among the independent factors. 


8 



The plots of studentized residuals and actual buckling load responses versus the predicted 
buckling loads for the fully quadratic RS model are shown in Fig. 3. 


Table 3 Comparison of candidate models for GRST 


Sample Size 

Sample Mean 

RMSE 

R 2 

Adj. R 2 

COV(%) PRESS/ SS Error 

300 

145372 

4684.0 

0.998 

0.976 

3.22 

194.0 

500 

143275 

4379.0 

0.991 

0.979 

3.06 

5.6 

1000 

143089 

4794.8 

0.981 

0.974 

3.35 

2.0 

1500 

143019 

4907.9 

0.978 

0.973 

3.43 

1.6 

2000 

143413 

4895.1 

0.977 

0.974 

3.41 

1.4 

2500 

143337 

4967.2 

0.976 

0.973 

3.47 

1.3 

3000 

143030 

4939.7 

0.976 

0.973 

3.45 

1.2 

3000L 

143030 

9464.0 

0.902 

0.901 

6.62 

1.02 

3000LQ 

143030 

8897.5 

0.914 

0.912 

6.22 

1.03 


The search for influential observations through the use of Cook’s D and Hat Diagonal did not 
indicate the presence of any outliers. 

To complete the examination of candidate RS models, each model was tested for accuracy based 
on its predictability of buckling response at design points not included in the multiple linear 
regression modeling. In this case, a set of 50 design points was generated using Latin Hypercube 
Sampling (LHS) feature of VisualDOC (2002) with the exact and predicted buckling loads 
evaluated at each design point. Table 4 gives a summary of results found in this examination. 
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(a) 



(b) 


Fig. 3 Plots of studentized residual (a) and actual buckling load (b) 
for the GRSM based on 3000 response samples 

These results indicate no statistical difference between the exact mean and the approximate mean 
predicted by each candidate GRS model. Although, in terms of maximum residual, the 3000- 
sample model performs the best, there is no significant difference among the models beyond 
1000 samples; therefore, any of these models would be an acceptable candidate for GRST. Upon 
considering all the available information, we opted to use the model based on 3000 response 
samples for GRST-based solution of the RBO problem. 
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T able 4 Predictability of candidate GRS m odels 


Model 

Per, lb 

^max> Ik 

Eave. lb 

Exact 

138535 

— 

— 

500 

139684 

22569 

6713 

1000 

139163 

12361 

4228 

1500 

138260 

12001 

4537 

2000 

139984 

12357 

4652 

2500 

139902 

12225 

4738 

3000 

139423 

11739 

4651 


A computer program based on the procedure described in Fig. 2 was developed to obtain the 
solutions for the RBO problem. The results presented in Table 5 below are based on the 
following conditions: Mean ply angles (0j / ... / d$) s = (±45/ 90 4 /+45) s ; Ply angle standard 

deviation Oq, = 3 deg.; Side constraints y} =0.0029 in., Y“ =0.0065 in.; Initial values of design 
variables, Yj =0.005 in.; Three combinations of mean length and diameter; Two different values 
for COV of the applied load; Target reliability: fi min = 3.72 (i.e., P f = 0.0001) 


The quantities shown are the mean values for the estimated buckling load, P cr , cylinder weight, 
W, and cylinder wall thickness, h. 


Table 5 RBO results 

using GRST 


Mean value at optimum design 

Parameter 

COV(P fl ) = 5% 

' COW(P a ) = 10% 

Mean diameter = Mean length =10 in. 

Per , lb 

124,537 

141,517 

W. lb 

1.31 

1.40 

h, in 

0.0724 

0.0778 

Mean diameter = 10 in.. Mean length = 20 in. 

Per , lb 

128,560 

143,969 

W, lb 

2.73 

2.92 

h, in 

0.0757 

0.0808 

Mean diameter = 20 in., Mean length = 10 in. 

Per, lb 

126,016 

142,396 

W, lb 

2.57 

2.75 

h, in 

0.0715 

0.0765 


The effect of uncertainty in the applied load is illustrated through the use of two different COV 
values (i.e., 5% and 10%). 

3.2 Sequential local response surface technique (SLRST) 

Although GRST provides considerable computational time saving by eliminating repeated exact 
buckling analyses in the evaluation of (3 b and the solution of RBO problem, it still requires a 
rather large number of simulation cycles for development of an accurate GRSM at pre- 
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optimization stage. A potentially more efficient alternative to GRST is the sequential local 
response surface technique (SLRST). 

In SLRST, the design space is partitioned into smaller subregions and in each subregion a linear 
RS model (LRSM) in the form 

P cr (X) = a 0 + la i X i (10) 

i-l 


is used to approximate the axial buckling force. 

Inspired by the sequential local approximate optimization technique, the procedure for SLRST, 
shown in Fig. 4, is similar in principle to the multipoint approximation method of Toropov et al. 
(1993) and Polykin et al. (1995), and the sequential local optimization strategy implemented by 
Sevant et al. (2000). 



Fig. 4 Flow chart of the optimization procedure using SLRST 

The main advantage of SLRST is that since the limit state function is described by a linear 
equation, the reliability index can be calculated using the basic reliability formula 


Pb = 




( 11 ) 


where up , fXp a are the mean values of P cr and P a , respectively, and op , a/> are the 
corresponding standard deviations with the variance of P cr determined as 
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( 12 ) 



where the partial derivatives of P cr are found analytically from the LRSM. 

Because in SLRST, a separate LRSM is developed for each subregion optimization, the overall 
efficiency of this approach is governed by how efficiently the unknown coefficients of each 
LRSM are calculated. 

With 22 random variables, Eq. (10) has a total of 23 unknown coefficients. Theoretically, only 
23 design points and corresponding response values are needed for the calculation of the 
unknown coefficients. However, is it possible to find a reasonably accurate LRSM with only 23 
response observations? Furthermore, does the accuracy of LRSM depend on how the design 
points are selected? 

To answer these questions, we conducted a preliminary study by developing six separate LRS 
models based on six different experimental and random design approaches (i.e., Koshal, Box- 
Behnken, Taguchi, Latin Hypercube Sampling, Random, and Quadrature). 

The Random design points are obtained using a random number generator while the Quadrature 
design points are found using the non-uniformly spaced quadrature formulae described in 
Appendix B. Overall, Taguchi and Box-Behnken require 81 and 44 design points, respectively, 
while the remaining four approaches require only 23. 

Table 6 shows the listing of random variables, their mean values as well as the corresponding 
lower and upper bounds used in developing the six candidate LRS models. The lower and upper 
bounds in Table 6 are found using the quadrature formulae in Appendix B and could be 
interpreted as those defining a single subregion in SLRST implementation. 

With the exception of Quadrature design, the design points for other approaches were obtained 
using VisualDOC (2002). Furthermore, the unknown coefficients in all LRS models were 
obtained using the REG procedure of SAS. 

To examine their accuracy, we used each candidate LRSM, defined in Table 7, to calculate the 
buckling load with all random variables kept at their respective mean values. These results are 
shown in column 2 of Table 8. In addition, we generated 25 random design points (other than 
those used in developing any of the LRS models but consistent with the bounds in Table 6) and 
calculated the response value at each location using the six LRS models. Columns 3 and 4 in 
Table 7 show the maximum and average residuals for each model with the ratio of average error 
to exact mean buckling load shown in column 5. 

It is clear from these results that the selection of design points does make a difference in the 
accuracy of the resulting LRSM as Taguchi, Box-Behnken, and Quadrature designs produced 
more accurate models than the rest. However, in terms of efficiency, the Quadrature design is 
the most superior, as it requires only 23 design experiments compared to 81 for Taguchi and 44 
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for Box-Behnken. Therefore, based on these findings, the Quadrature design was subsequently 
adopted for the SLRST-based solution to the RBO problem. 

Table 6 Bo unds on random variables according to the quadrature formulae in A ppendix B 


Random Variable Mean Value Lower Bound Upper Bound 


£,,psi 

1.800E+07 

1.795E+07 

1.912E+07 

E 2 , psi 

1.350E+06 

1.345E+6 

1.456E+06 

v I2 

0.226 

0.225 

0.247 

^12> 

5.430E+05 

5.403E+05 

5.939E+05 

t\, in 

0.005 

0.0049 

0.0061 

t 2 , in 

0.005 

0.0049 

0.0061 

h, in 

0.005 

0.0049 

0.0061 

t 4 , in 

0.005 

0.0049 

0.0061 

*5. in 

0.005 

0.0049 

0.0061 

h, in 

0.005 

0.0049 

0.0061 

f 7 , in 

0.005 

0.0049 

0.0061 

h, in 

0.005 

0.0049 

0.0061 

0„ deg 

45.0 

40.8 

58.7 

0 2 » deg 

-45.0 

-48.2 

-31.4 

0 3 , deg 

90.0 

87.4 

103.6 

0 4 * deg 

90.0 

87.4 

103.6 

0s, deg 

90.0 

87.4 

103.6 

0 6 . deg 

90.0 

87.4 

103.6 

0 7 . deg 

-45.0 

-48.2 

-31.4 

0 g , deg 

45.0 

40.8 

58.7 

D, in 

20.0 

18.0 

24.0 

L, in 

10.0 

8.3 

11.7 


An example of the sequence of subregion optimization is depicted in Fig. 5 with the initial 
design pointed designated by S. With COV of each ply thickness held fixed, the size of each 
subregion is determined based on the subregional (local) lower and upper bounds on the design 
variables defined as 


vhU 

*i,k 


Yi,k *6% 


i = 1,2,..., 8 


(13) 


where is the value of ith design variable at the beginning of kth subregion optimization cycle 
with the corresponding lower and upper move limits calculated as 



(14-a) 


(14-b) 
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Table 7 Coefficients in each candidate LRSM 


Random 

Variable 

RS 

Coefficient 

Koshal 

B-B 

Taguchi 

LHS 

Random 

Quadrature 

— 

Oo 

-62107.69 

-191034.77 

-188337.22 

806970.26 

-422810.24 

-250237.00 

£, 


0.006402 

0.008339 

0.006432 

-0.066344 

0.014638 

0.008319 

E 2 

a 2 

0.011382 

0.015254 

0.006039 

0.143120 

-0.007214 

0.026950 

V l2 

(X^ 

5909.09 

6818.18 

-319.87 

6517713.61 

-156591.33 

85139.52 

On 

a 4 

0.024231 

0.0328 

0.021983 

-2.2865 

0.012186 

0.058529 

t\ 

a 5 

5730769 

7446154 

8329630 

12846767 

13073274 

7366942 

h 


9800000 

11161538 

10094302 

-11496414 

11840124 

10705970 

h 


7276923 

9069231 

7733333 

-16173878 

1542889 

8878919 

u 


7276923 

9069231 

7118234 

7367939 

10791940 

8987386 

h 

a g 

7276923 

9069231 

7927350 

-18894534 

-67322 

9107382 

1 6 

a \o 

7276923 

9069231 

8138177 

26805393 

-10111602 

9246309 

h 

flu 

5476923 

6000000 

7515385 

74168144 

6383828 

7808193 

h 

a i2 

6161538 

6753846 

8390313 

36139560 

16876669 

8444851 

0. 

fll3 

-1119.55 

-1086.59 

-1148.46 

2111.13 

-323.175 

-1093.78 

02 


1032.14 

1706.55 

1320.66 

-1048.38 

872.72 

1034.13 

03 

flis 

-350.000 

-239.506 

-68.930 

-1302.09 

455.058 

-280.605 

04 

fll6 

-387.037 

-251.235 

-121.491 

1033.4 

-66.6697 

-301.852 

05 

a \l 

-417.901 

-259.877 

-122.634 

1 106.4 

288.169 

-315.524 

06 

a \% 

-441.975 

-268.519 

-159.694 

-4298.37 

604.96 

-320.028 

07 

a \9 

223.214 

167.857 

14.3078 

-3174.94 

587.593 

68.8512 

08 

a 20 

337.989 

482.682 

206.683 

1031.57 

1399.56 

162.724 

D 

a 2\ 

18.3333 

-608.333 

105.062 

-20840.49 

-1387.45 

-318.352 

L 

®21 

2029.41 

2802.94 

1572.88 

-123.604 

1193.18 

2069.58 


The 8 thickness design variables correspond to the mean values of random variables 5 through 
12. In that group, the twelfth random variable has the smallest perturbations away from the 
mean according to Quadrature design (see Appendix B). Therefore, the move limits in Eq. (14) 
are specified based on perturbations in the twelfth random variable. 


Table 8 Comparison of candidate LRS models 



fip„ lb e m ax> lb 

^ave’ ^ 

£ avJ Pcr exaa 

Exact 

1.543e5 - 

— 

— 

Koshal 

1.497e5 2.17e4 

8.58e3 

4.6% 

B-B 

1.498e5 1.72e4 

5.34e3 

2.9% 

Taguchi 

1.465e5 1.21e4 

4.27e3 

2.3% 

LHS 

1.932e5 1.76e5 

8.25e4 

44.4% 

Random 

1.361e5 2.93e4 

1.27e4 

6.8% 

Quadrature 1.528e5 1.22e4 

6.00e3 

3.2% 


For COV(f) = 5% and n = 22, we find 6 l iJc = 0.02087 Y iJk and d“ k = 0.22958 Y i k , which indicate 
unsymmetric move limits with more than a 10 to 1 ratio favoring freedom of movement towards 
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the upper bound. If at the end of &th subregion optimization cycle, design variables have 
different values, then they will be subject to different move limits in cycle k+l. 



Fig. 5 A notional schematic of SLRST over four subregions for a 2-dimensional design space 

By staying within the specified move limits, the buckling load estimates obtained from LRSM 
are found to be fairly accurate with an observed error of no more than roughly 1.5 to 2%. 

With the local bounds defined, design improvement is sought in the current subregion by solving 
the mathematical programming problem in Eq. (6). Upon completion of one subregion 
optimization cycle, the search continues to the next cycle focusing on an updated subregion with 
location defined by the final values of design variables in the current subregion and its size 
determined by Eqs. (13) and (14). 

The move limits are reduced in half if the absolute difference between objective functions in two 
consecutive cycles is less than 6x1 O' 4 , and hard convergence is reached when the difference is 
less than or equal to 6x1 O' 6 . 

The SLRST-based optimization results for the same cylinder geometry, loading conditions, and 
target reliability as in Table 5 are shown in Table 9. 

These results are based on an initial mean ply thickness (design variable) of 0.004 in. From this 
initial design, it takes an average of 10 subregion optimization cycles to reach optimal design. 
Because of unsymmetric move limits, the choice of initial design can influence the number of 
subregion optimization cycles performed. To demonstrate this characteristic, the convergence 
history for E[D]/E[L]/COV(PJ = 20/10/10 cylinder with two different initial design points is 
shown in Fig. 6. Curve A is for the initial mean ply thickness of 0.006 in. while Curve B is for 
0.004 in. 
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Table 9 RBO results using SLRST 



Mean value at optimum design 

Parameter 

CO V(P a ) = 5% COV(PJ = 10% 

Mean diameter = Mean length =10 in. 

P cr , lb 

129463 142953 

W. lb 
h, in 

1.33 1.42 

0.0737 0.0789 


Mean diameter = 10 in., Mean length = 20 in. 


lb 

129375 

145163 

W, lb 

2.78 

2.96 

h, in 

0.0771 

0.0821 

Mean diameter 

= 20 in., Mean length = 10 in. 

Per, lb 

132311 

147664 

W, lb 

2.63 

2.79 

h, in 

0.0732 

0.0776 


Since the thickness in each ply is treated as an independent random variable, it is possible for the 
plies that have the same mean orientation angles to have different thickness values. Figure 7 
shows the variation in the eight design variables from initial to final point corresponding to 
Curve B in Fig. 6. 



Fig. 6 Convergence history based on two different initial design points 

To validate the accuracy of our reliability analysis and to determine whether the optimal design 
meets the specified target reliability index and probability of failure, we performed an 
independent analysis using the structural reliability code, NESSUS (2002). Because of its 
reasonably good efficiency and accuracy, we chose the analytical method AMV+ for reliability 
index and probability of failure calculations. In this case, NESSUS calls the exact buckling 
analysis code directly for each limit state function evaluation and has no interaction with any of 
the RS models developed here. Table 10 shows the reliability index and probability of failure 
values for the optimal designs obtained through SLRST. 
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Fig. 7 Convergence history of design variables for the cylinder with 
E[D]/E[L]/COV(P a ) = 20/10/10 


To measure the effect of ply pattern on optimal design, SLRST was used to obtain optimal 
design solutions for cylinders with four different mean ply orientation angles. The results shown 
in Table 1 1 are for the E[D]/E[L]/COV(FJ = 20/10/10 cylinder with target reliability index of 
3.72 and mean ply orientation angles and stacking sequence as follows: PP1: (±45/ 9 O 4 /+45) s ; 

PP2: (±45 / 0 4 / +45 ) s ; PP3: (±45 / 0 / 90) 2s ; PP4: (0 / ±30 / ±45 / ±60 / 90 ) s 


Parameter 

SLRST 

NESSUS 

Mean diameter = Mean length =10 in. 

P 

3.72 

3.78 

P f 

1.0e-4 

0.78e-4 

Mean diameter = 

10 in.. Mean length = 20 in. 

P 

3.72 

3.63 

P f 

1.0e-4 

1.4e-4 

Mean diameter = 

20 in.. Mean length = 10 in. 

P 

3.72 

4.16 

P f 

1.0e-4 

0.16e-4 


Table 1 1 SLRS T-based RBQ results for cylinders with various me an ply patterns 


Parameter 

PP1 

PP2 

PP3 

PP4 

P^.lb 

147664 

146653 

146320 

148682 

W, lb 

2.79 

2.87 

2.59 

2.56 

h, in 

0.0776 

0.0799 

0.0720 

0.0712 


Figure 8 shows the convergence history for the designs in Table 1 1 with the optimal mean ply 
thickness values shown in Table 12. 
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Optimization Cycle 

Fig. 8 Convergence history for cylinders with four different ply patterns 


Table 12 SLRST-based RBQ results for cylinders with various mean pl y patterns 


Mean Ply Thickness PP1 PP2 PP3 PP4 



in 

4.4E-03 

3.3E-03 

3.3E-03 

3.6E-03 

tl* 

in 

6.2E-03 

3.3E-03 

4.3E-03 

3.6E-03 


in 

4.6E-03 

4.7E-03 

4.0E-03 

4.3E-03 

^ 4 > 

in 

4.9E-03 

4.9E-03 

4.4E-03 

5.3E-03 


in 

5.2E-03 

5.2E-03 

4.9E-03 

4.4E-03 

h. 

in 

5.5E-03 

5.4E-03 

5.2E-03 

5.0E-03 

* 7 . 

in 

3.9E-03 

6.5E-03 

5.0E-03 

4.5E-03 


in 

4.2E-03 

6.5E-03 

5.1E-03 

4.9E-03 


The SLRST and GRST-based solutions are comparable in terms of accuracy. However, in terms 
of computational efficiency, the two techniques are considerably different. While in GRST, the 
number of exact buckling response evaluations is independent of the number of optimization 
iterations, in SLRST it is not, as each subregion optimization cycle requires 23 exact buckling 
load evaluations. Depending on the model used, GRST required from 1000 to 3000 design points 
and exact buckling response evaluations as indicated in Tables 3 and 4. By contrast, the number 
of design points and responses required by SLRST was 23 times the number of optimization 
cycles, which varied from a minimum of 3 to a maximum of 20 for a total of 69 to 406, 
respectively. 


4. Reliability Sensitivity Analysis 

The probabilistic sensitivity derivatives of reliability index in Eq. (5) are evaluated at the design 
point x * and normalized as 





i = l,2,...,n 


(13-a) 
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i = 12 


(13-b) 


Vi - 


dh 


do 


X> 



and plotted in Fig. 9 for a composite cylinder with (±45/ 9 O 4 /+45) s mean ply pattern and 
E[Z3]=E[L]=15 in. In this case, the GRSM is used for the evaluation of the sensitivities. 





0.0 T 


-0.2 

.0.4 

-0.6 

-0.8 

-1.0 -* 

(b) 

Fig. 9 Normalized sensitivity derivatives of fa with respect to the mean (a) 
and standard deviation (b) of each random variable 

Of all random variables, the applied load is found to have the greatest influence on fa followed 
by elastic modulus in fiber direction and cylinder diameter. 

Although total shell thickness has a large impact onP cr , the influence of each individual ply 
thickness is less significant. 

As for the effect of parametric uncertainty, it is evident from Fig. 9-b that the amount of scatter 
in the applied load has the largest effect on fa. Without the sensitivity derivatives of fa with 
respect to P a , Fig. 9 is replaced by Fig. 10, which more clearly shows the influence of random 
variables on fa through P cr . 
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max 


1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

- 0.2 



(a) 



Fig. 10 Normalized sensitivity derivatives of with respect to the mean (a) 
and standard deviation (b) of random variables affecting Per 

5. Concluding remarks 

This report presented the development of global and sequential local response surface techniques 
for reliability-based optimization of composite structures. These techniques were subsequently 
applied to reliability-based optimization of a composite circular cylinder under axial 
compression and buckling instability. 

Based on the formulation and solution of the cylinder RBO example problem, the following 
conclusions are reached: 

• As in most RBO problems, the computational requirement for reliability index evaluation 
far exceeds that for design space search of optimal point. 

• RBO results using GRST and SLRST have the same level of accuracy. 

• SLRST has a superior computational efficiency than GRST as in most cases it required 5 
to 10 times fewer exact buckling response evaluations. 

• The choice of experimental design technique affects the accuracy and efficiency of LRS 
models with the Quadrature design found to be most efficient. 

• The non-uniformly spaced quadrature formulae result in greater freedom of movement 
toward the upper than the lower bound in each subregion optimization. This could speed 
up or slow down the search toward optimal design depending upon the choice of initial 
design point. 

• Among all random variables, the applied load followed by elastic modulus in the fiber 
direction and diameter were found to have the greatest influence on cylinder reliability 
index. 
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• Besides the mean value, the amount of scatter in applied load is found to have a large 
influence on cylinder reliability and optimal cylinder wall thickness. 

• The choice of laminate ply pattern affects cylinder reliability through changes in the 
buckling load. 
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Appendix A 
Shell buckling analysis 

The computational model used in buckling analysis of circular cylinders is a semicircular shell 
with the unloaded edges having symmetric boundary conditions and the loaded edges treated as 
clamped as shown in Fig. Al. This model falls under the general category of thin-shell 
structures with the displacement field described by the first-order shear deformation theory 
formulated as 

u(x,y,z) = u 0 (x,y) + z<p x (x,y ) 

v(x,y,z)= v 0 (x,y)+ z(p y (x,y) (Al) 

w(x,y,z) = WoC^y) 

where Ho^O’^O are the midplane displacements in x, y, z directions (see Fig. Al), respectively, 
and <(> x ,(py describe rotations about the y and x axes, respectively. 

The strain-displacement relations are based on Sanders-Koiter shell theory with the in-plane 
strains (e*,^,/*^,, at z = 0), transverse shear strains (Yxz'Yyz* at z = 0), and curvatures 
{k x ,k y ,K xy) formulated as 
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dx 
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dy 
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dx 
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d( Px | , W ^0 , 

dx 2/?\ dy dx) 


0 , 

Yxz =< Px + 


dw Q 
dx 


Y 0 = <* + ^0- 


10 

/? 


(A2) 


where R is the shell radius of curvature in y direction (see Fig. Al). Hence, the elastic strain 
energy stored in the shell can be expressed as 


u-\Sa{°} t 


Ay A/ ® 


sym 


Du 0 


-pq 


{ e}dA 


(A3) 


where A i j,Bij,D i j are the laminate extensional, coupling, and bending stiffness matrices, 
respectively, and C pq is the transverse shear stiffness matrix with the strain vector defined as 

{e} r = {e? ej y °xy k x K y Kxy Yxz y (A4) 

The work done by the applied edge load (see Fig. Al) is calculated using the nonlinear 
components of strain-displacement relations with the resulting equation expressed as 


w - \s a n x 


(*S£\ 2 

\ dx ) 


JdWQ^' 
V dx ) 


dA 


(A5) 


where N x is the stress resultant in the x direction. 
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To obtain the critical buckling load, the displacements mq^O’^O anc ^ rotations (f> x ,(py are 

approximated by different Ritz series with Legendre polynomials used as the interpolation 
functions such that the essential boundary conditions are satisfied. Then through the application 
of the principle of minimum total potential energy the critical buckling load is found by setting 
Eqs. (A3) and (A5) equal to each other and solving the resulting eigenvalue problem for the 
critical load factor such that 

N cr = X cr N x (A6) 






Fig. A1 Computational model of circular cylinder 

A computer implementation of the described analysis procedure developed by Jaunky and 
Knight (1999) is used to generate the necessary response samples for development of GRS and 
LRS models. 

Prior to its use, we compared the results obtained from this code with those from linear buckling 
analysis in MSC/NASTRAN. The latter solutions were based on a 60 x 216 finite element 
model of the complete cylinder with approximately 76,700 total degrees of freedom. The critical 
buckling loads obtained from these two programs were found to be in excellent agreement. 

Appendix B 

Selection of design points based on quadrature formulae 

The approach we refer to as Quadrature design for fitting a linear RS model is based on targeted 
sampling of design space at n + 1 points corresponding to n+1 vertices of a regular n-simplex with 
centroid at the origin. The coordinates of ith design point are obtained using the relationship 

*i,j = PXj + 3 J °X } * - 1,2, ...,n +1 (Bl) 

j =1,2 ,...,n 

where / Ux.,°Xj represent the mean value and standard deviation of y'th random variable, 
respectively, and Zjj is the perturbation vector at ith design point. 
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The elements of Zjj correspond to the quadrature points developed by Stroud and Secrest (1963) 
for approximate solution of multiple integrals of polynomial functions involving two or more 
independent variables, and later used by Zhou and Nowak (1998) for non-product point 
integration of a function of multiple normal random variables. 


At the first design point, the perturbation vector is given as 

Zij -(ZU’^’-ZU) = (VM>,0,...,0) (B2) 

Therefore, while the first design variable is perturbed from its mean value by Jnox , , the 

remaining variables are kept at their respective mean values. For the subsequent n design points 
(see Fig. Bl), the perturbation vector is defined as 


„) 

•) 


7 ■=(- II IZH I (n + 0 I ( ” + l) 

,J { \n’ | n(n - 1) ’ \ (n - l)(n - 2) ’ \ (n - 2)(n - 3) ’ 


z. ( II [FIE I ( w+1 >~ I (» + o~ 

^ l+1J ( Vn’ \ n(n - 1) ’ y(n-l)(n-2)’ \ (n - 2)(n - 3)" 



(B3) 



Fig. Bl Depiction of random values selected for design point 
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